R Notebook
L’objectif de ce TP est de traiter le cas d’un mélange gaussien (à deux composantes pour le moment).
Simulation d’un mélange
On peut commencer par simuler un mélange pour travailler avec des données “faciles” :
rmix = function(N,p,moy1,moy2,t1,t2){
b = rbinom(N,1,p)
x = (1-b)*rnorm(N,moy1,1/sqrt(t1))+(b)*rnorm(N,moy2,1/sqrt(t2))
return(x)
}
hist(rmix(1000,0.6,0,5,1,1),100,main = "Simulated Mixture")
Modèle à variables cachés
On travaille maintenant avec le modèle suivant :
On peut initialiser les données :
N = 1000
true = list(p = 0.6,m1 = 0,m2 = 4,tau1 =1,tau2 = 5)##pour la simulation
x = rmix(1000,true$p,true$m1,true$m2,true$tau1,true$tau2)
h = list(a =1,b = 1, mu = 0,sig = 1000,k = 1,lam = 0.001)
l = list(p = 1/2,m = c(0,0),tau = c(1,1),z = 1+rbinom(length(x),1,1/2))init = function(x){
h = list(a =1,b = 1, mu = 0,sig = 1000,k = 1,lam = 0.001)
l = list(p = 1/2,m = c(0,0),tau = c(1,1),z = 1+rbinom(length(x),1,1/2))
return(list(h = h,l = l))
}
logvraiss= function(l,x){
return( sum( log( dnorm(x,l$m[l$z],1/sqrt(l$tau[l$z]) ) ) ) )
}
logvraiss(l,x)[1] -5851.625logprior = function(l,h){
val = sum(log(dbinom(l$z-1,1,l$p )))+
log(dbeta(l$p,h$a,h$b))+
sum(log(dnorm(l$m,h$mu,h$sig) ))+
sum(log(dgamma(l$tau,h$k,h$lam)))
return(val)
}
logprior(l,h)[1] -722.6181logpost = function(l,x,h){
return(logvraiss(l,x)+logprior(l,h))
}
logpost(l,x,h)[1] -6574.243step_MH = function(l,x,h,sigprop = c(0.01,0.01,0.01),nrunz = 50){
##variance prop codées en dur pour le moment
sig_prop_m = sigprop[1]
inv_prop_tau = 1/sigprop[2]
inv_prop_p = 1/sigprop[3]
lnew = l
N = length(x)
# print(logpost(lnew,x,h))
##ici, c'est bien couteux de recalculer tout à chaque fois -> on devrait calculer que pour x[i] la vraisemblance !!!!
##proposition m
i = sample(2,1)
lnew$m[i] = lnew$m[i]+rnorm(1,0,sig_prop_m)
if (log(runif(1))>(logpost(lnew,x,h)-logpost(l,x,h)) ){
##on rejette avec proba 1-alpha, ici la proposition est symétrique
lnew$m[i] = l$m[i]
}
# print(logpost(lnew,x,h))
if (lnew$m[2]<lnew$m[1]){
##si les m se sont croisés, on échange tout : m, tau, p et z
lnew$m = rev(lnew$m)
lnew$tau = rev(lnew$tau)
lnew$p = 1-lnew$p
lnew$z = 3-lnew$z
}
# print(logpost(lnew,x,h))
##proposition tau
i = sample(2,1)
lnew$tau[i] = rgamma(1,inv_prop_tau*l$tau[i]**2,inv_prop_tau *l$tau[i])
log_ratio_prop = log(dgamma(l$tau[i] ,inv_prop_tau*lnew$tau[i]**2,inv_prop_tau *lnew$tau[i]))-log(dgamma(lnew$tau[i] ,inv_prop_tau*l$tau[i]**2,inv_prop_tau *l$tau[i]))
if (log(runif(1))>(log_ratio_prop+logpost(lnew,x,h)-logpost(l,x,h)) ){
##on rejette avec proba 1-alpha. ici on prend en compte le ratio
lnew$tau[i] = l$tau[i]
}
##proposition z
# for (i in 1:N){
# nrunz = 50
samp = sample(N,nrunz)
for (i in samp){
# if (runif(1)<1/2){
# i = sample(N,1)
lnew$z[i] = 3-l$z[i]
if (log(runif(1))>(logpost(lnew,x,h)-logpost(l,x,h)) ){
##on rejette avec proba 1-alpha. ici la proposition est symétrique
lnew$z[i] = l$z[i]
} else{
l$z[i] = lnew$z[i]
}
# }
}
# print(logpost(lnew,x,h))
##proposition p
lnew$p = rbeta(1,inv_prop_p*l$p,inv_prop_p*(1-l$p))
log_ratio_prop = log(dbeta(l$p,inv_prop_p*lnew$p,inv_prop_p*(1-lnew$p)))-log(dbeta(lnew$p,inv_prop_p*l$p,inv_prop_p*(1-l$p)))
if (log(runif(1))>(log_ratio_prop+logpost(lnew,x,h)-logpost(l,x,h)) ){
##on rejette avec proba 1-alpha. ici on prend en compte le ratio
lnew$p = l$p
}
return(lnew)
}
lnew = step_MH(l,x,h)run_MH = function(nMH,x,...){
linit = init(x)
l = linit$l
h = linit$h
traj = matrix(0,nMH,5)
colnames(traj) = c("p","m0","m1","tau0","tau1")
nb = rep(0,length(x))
for (i in 1:nMH){
nb = nb+(l$z-1)
l = step_MH(l,x,h,...)
traj[i,] = c(l$p,l$m,l$tau)
}
class = nb/nMH
return(list(traj = traj ,class =class,l = l))
}
nMH = 100
ltraj = run_MH(nMH,x)
plot(1:nMH,ltraj$traj[,2],"l")
lines(1:nMH,ltraj$traj[,3],col = "red")
Pour un plot plus joli, on reprend la fonction de la dernière fois
plot_run = function(x,mem,l,i = NULL,nbpoints = 10000){
if (is.null(i)){i = dim(mem)[1]}
par(mfrow = c(1,2))
hist(x,50,freq= F,main = "Current estimated density")
curve((1-l$p)*dnorm(x,l$m1,1/sqrt(l$tau1))+(l$p)*dnorm(x,l$m2,1/sqrt(l$tau2)),add = T,col = "red",lty = 2,lwd= 2)
curve((1-l$p)*dnorm(x,l$m1,1/sqrt(l$tau1)),add = T,col = "blue",lty = 3,lwd = 2)
curve((l$p)*dnorm(x,l$m2,1/sqrt(l$tau2)),add = T,col = "green",lty = 3,lwd = 2)
z1 = sample(i,nbpoints,replace = T)
z2 = sample(i,nbpoints,replace = T)
plot(z1,mem[z1,2]+2*1/sqrt(mem[z1,4])*runif(nbpoints,-1,1),ylim = range(x),col = "cyan",cex = 0.1,main = "MH trajectory with 2-sigma bands",ylab = "Trajectories for m1,m2",xlab= "Time")
points(z2,mem[z2,3]+2*1/sqrt(mem[z2,5])*runif(nbpoints,-1,1),ylim = range(x),col = "yellow",cex = 0.1)
lines(1:i,mem[1:i,2],ylim = range(x),col = "blue")
lines(1:i,mem[1:i,3],col = "green")
}
init = function(x){
h = list(a =1,b = 1, mu = 0,sig = 1000,k = 1,lam = 0.001)
l = list(p = 1/2,m = c(0,0),tau = c(1,1),z = 1+rbinom(length(x),1,1/2))
# l = list(p = 1/2,m = c(-1,1),tau = c(1,1),z = 1+(x>mean(x)))
return(list(h = h,l = l))
}
N = 1000
true = list(p = 0.6,m1 = 0,m2 = 2,tau1 =1,tau2 = 5)##pour la simulation
x = rmix(1000,true$p,true$m1,true$m2,true$tau1,true$tau2)
nMH = 2000
sigprop = c(0.5,0.1,0.01)
nrunz = 50
ltraj = run_MH(nMH,x,sigprop,nrunz)
l = ltraj$l
lend = list(p = l$p,m1 = l$m[1],m2 = l$m[2],tau1 = l$tau[1],tau2 = l$tau[2])
plot_run(x,ltraj$traj,lend)
On peut alors voir les paramètres finaux :
cat("estimés : ",unlist(lend),"\n")estimés : 0.5321751 0.4869726 1.986869 0.7883785 5.634391 cat("théorique :",unlist(true),"\n")théorique : 0.6 0 2 1 5 On peut aussi observer la propotion du temps que chaque observation a passé avec :
s = sort(x,index.return = T)
plot(s$x,ltraj$class[s$ix],"l",main = "Classification probability",xlab = "x",ylab = "Estimated P(Z(x)= 1)")
Et si on veut faire un petit film, il n’y a qu’à modifier les fonctions
library(animation)
plot_run = function(x,mem,class,l,i = NULL,nbpoints = 10000,MHtot = NULL){
if (is.null(i)){i = dim(mem)[1]}
if (is.null(MHtot)){MHtot = dim(mem)[1]}
# par(mfrow = c(1,3))
hist(x,50,freq= F,main = "Current estimated density")
curve((1-l$p)*dnorm(x,l$m1,1/sqrt(l$tau1))+(l$p)*dnorm(x,l$m2,1/sqrt(l$tau2)),add = T,col = "red",lty = 2,lwd= 3)
curve((1-l$p)*dnorm(x,l$m1,1/sqrt(l$tau1)),add = T,col = "blue",lty = 3,lwd = 3)
curve((l$p)*dnorm(x,l$m2,1/sqrt(l$tau2)),add = T,col = "green",lty = 3,lwd = 3)
z1 = sample(i,nbpoints,replace = T)
z2 = sample(i,nbpoints,replace = T)
plot(z1,mem[z1,2]+2*1/sqrt(mem[z1,4])*runif(nbpoints,-1,1),ylim = range(x),xlim = c(1,MHtot),col = "cyan",cex = 0.1,main = "MH trajectory with 2-sigma bands")
points(z2,mem[z2,3]+2*1/sqrt(mem[z2,5])*runif(nbpoints,-1,1),ylim = range(x),col = "yellow",cex = 0.1)
lines(1:i,mem[1:i,2],ylim = range(x),col = "blue")
lines(1:i,mem[1:i,3],col = "green")
s = sort(x,index.return = T)
plot(s$x,class[s$ix]/i,"l")
}
run_MH = function(nMH,x,pas_anim,...){
linit = init(x)
l = linit$l
h = linit$h
traj = matrix(0,nMH,5)
colnames(traj) = c("p","m0","m1","tau0","tau1")
nb = rep(0,length(x))
traj[1,] = c(l$p,l$m,l$tau)
for (i in 2:nMH){
nb = nb+(l$z-1)
l = step_MH(l,x,h,...)
lnow = list(p = l$p,m1 = l$m[1],m2 = l$m[2],tau1 = l$tau[1],tau2 = l$tau[2])
if ((i%%pas_anim)==2){
plot_run(x,traj[1:i,],class = nb,l = lnow,i=i,MHtot =nMH )
ani.record()
}
traj[i,] = c(l$p,l$m,l$tau)
}
class = nb/nMH
return(list(traj = traj ,class =class,l = l))
}
#
# N = 100
# true = list(p = 0.6,m1 = 0,m2 = 2,tau1 =1,tau2 = 5)##pour la simulation
# x = rmix(1000,true$p,true$m1,true$m2,true$tau1,true$tau2)
# sigprop = 0.1*c(0.5,0.1,0.01)
# nrunz = 20
# nMH = 500
# pas_anim = 10
# par(mfrow = c(1,3))
# ani.options(interval = 0.05, nmax = floor(nMH/pas_anim))
# ani.record(reset = TRUE)
# ltraj = run_MH(nMH,x,pas_anim= pas_anim,sigprop,nrunz)
# saveGIF(ani.replay(),movie.name = "animation_MH.gif")